Multiscale simulations in simple metals: a density-functional based methodology 
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We present a formalism for coupling a density functional theory-based quantum simulation to a 
classical simulation for the treatment of simple metallic systems. The formalism is applicable to 
multiscale simulations in which the part of the system requiring quantum-mechanical treatment is 
spatially confined to a small region. Such situations often arise in physical systems where chemical 
interactions in a small region can affect the macroscopic mechanical properties of a metal. We 
describe how this coupled treatment can be accomplished efficiently, and we present a coupled 
simulation for a bulk aluminum system. 
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I. INTRODUCTION 

In performing computer simulations of complex phys- 
ical systems, a premium is placed on accuracy and effi- 
ciency Typically, one of these qualities can be improved 
at the expense of the other. In recent years, a new ap- 
proach has emerged that addresses a class of problems in 
which important small-length-scale phenomena are con- 
fined to a small region of the system but can have an 
impact on the behavior over a much larger scale. A typ- 
ical case is the tip of a crack, where localized chemical 
reactions may affect the strength of interatomic bonding, 
which in turn can influence in a dramatic way the macro- 
scopic mechanical properties of the solid. Such problems 
fall under the rubric of "multiscale" phenomena, requir- 
ing a treatment that addresses important aspects at each 
scale. The novel feature of this type of simulation is to 
use an accurate but computationally demanding method 
to treat the region of the system in which small-length- 
scale degrees of freedom are important, and a faster but 
less accurate method with the small-length-scale physics 
"coarse-grained" , to treat the rest of the system. 

Multiscale approaches rely on successfully coupling the 
two (or more) regions involved, which is referred to as 
seamless coupling. There is no single notion as to what 
constitutes a seamless coupling, but generally the cou- 
pling should be accomplished in such a way that the fic- 
titious boundary between the two regions, which only 
exists in the coupled simulation and not in the real sys- 
tem, does not introduce any physical consequences. For 
instance, recently several papers have dealt with the is- 
sue of ensuring that phonons are not reflected by the 
boundary between the two coupling methodsQ, Q- In 
the consummate multiscale method, the resulting ener- 
getics or dynamics is indistinguishable from what would 
result from a calculation with the accurate method ap- 
plied to the entire system. This ideal would be achieved 
only if the two simulation methods involved were seam- 
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lessly matched at the boundary, and further, only if the 
part of the system treated by the faster, less accurate 
method was indeed free of important small-scale physics. 
Another important but obvious characteristic of a good 
multiscale method is that the computational overhead of 
performing a coupled simulation is not significant. More 
specifically, the computation time for the coupled simula- 
tion should be on the order of computation time required 
for the accurate method to treat the small detailed re- 
gion, since the time required for the less accurate method 
to treat the rest of the system is typically several or- 
ders smaller; nothing is gained if the coupling is so costly 
that the coupled method takes as long as using the ac- 
curate method to treat the whole system. When the ap- 
proach requires coupling a quantum mechanical method 
to a classical method, additional complications arise be- 
cause of the presence of electronic degrees of freedom 
in the quantum mechanical region; thus boundary con- 
ditions on the electron wavefunctions must be imposed 
at the interface between the regions. Density functional 
theory (DFT) provides a significant simplification over 
more direct quantum mechanical methods in that the 
calculation of ground state energies and forces requires 
the minimization of a functional of the electron density 
p(r) only [3j ■ Thus, in principle boundary conditions need 
only be imposed on p(r). This statement only applies to 
the formulation of the problem that does not invoke the 
explicit calculation of electronic wavefunctions (the most 
common way of implementing DFT actually does involve 
individual electronic wavefunctions, the so-called Kohn- 
Sham orbitalsjj]). Coupling an approximate DFT calcu- 
lation that is based on the electronic density alone to a 
classical interatomic potential should be more straight- 
forward than coupling an orbital-based quantum mechan- 
ical method to a classical method. 

The present article describes a formalism for concur- 
rently coupling a system consisting of two regions, one 
treated with density functional theory (without invoking 
electronic wavefunctions) and the other with a classical 
interatomic potential. Due to the type of approximations 
involved, the present approach is particularly well suited 
for simple metallic systems; we emphasize, however, that 
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this is not an inherent limitation of our approach, but 
rather a limitation imposed by the shortcomings of the 
methodologies employed for the treatment of the vari- 
ous parts of the system, and if these are eliminated the 
approach could be generally applicable. In section [HI 
other methodologies for coupling multiple simulation ap- 
proaches and their relevance to the present methods are 
discussed. In Section IIIII the formalism of the present 
class of coupling methods is established. In Section llVl 
details of implementing the methods and achieving effi- 
ciency are presented and some tests of the method are 
reported in Section Finally we conclude in Section IVll 
with a discussion of the results of the tests. 



II. BACKGROUND 

A large number of concurrent multiscale methods 0,0, 
approach the problem of coupling two different simula- 
tion methods by writing the energy of the whole system 
as 

E[I + II] = E^I] + E 2 [II] + E int [1, 1 1] (1) 

where here and throughout the article, I refers to the 
small region where detailed physics are relevant, and II 
refers to the rest of the system. E±[I] represents the 
energy of region / with region II providing appropriate 
boundary conditions, E 2 [II] represents the energy of re- 
gion II in the same sense, and E[I + II] represents the 
total energy of the combined system. Eq. (JIJ expresses 
the total energy of the system as the energy of region I 
evaluated with the accurate simulation method 1, plus 
the energy of region II evaluated with the faster simu- 
lation method 2, plus an energy of interaction between 
the two subsystems. Typically, the crux of a multiscale 
method lies in its handling of E mt . Although tautologi- 
cal, Eq. (JTJ can be rearranged to yield an expression for 
the interaction energy: 

E int [1, 1 1] = E[I + II] - Ei[I] - E 2 [II] (2) 

This expression contains no new content, and merely 
serves to define E mt [I, II] , but nevertheless provides di- 
rection towards its calculation. 

The QM/MM methods are designed to achieve a goal 
similar to that of the present method, namely the cou- 
pling of a quantum-mechanical simulation with classical 
potentials, but in the context of covalently bonded or- 
ganic molecules. In such systems, bonds are localized 
and typically can be associated with two atoms at either 
end. The strategy often employed in QM/MM meth- 
ods to couple quantum mechanics to molecular mechan- 
ics is as follows the system is divided into QM and 
MM regions with a boundary that cuts across covalent 
bonds; Eqm is evaluated for the QM subsystem, plus 
additional 'link atoms' placed on the MM side of the sev- 
ered covalent bonds to mimic the system removed from 
the QM region; Emm is evaluated for the MM subsystem 



without the link atoms; and Eqm /mm consists of energy 
terms such as bond-bending terms that are left out of 
Eqm + Emm- A similar methodology was developed by 
Broughton et al.Q for quantum-classical coupling in sil- 
icon, the prototypical covalently-bonded bulk material. 
Such approaches rely on a somewhat artificial partition- 
ing of the total energy (e.g. into bond-bending and bond- 
stretching terms), and hence lack a definition that could 
be readily generalized. But due to the locality of physics 
in covalently-bonded systems for which QM/MM meth- 
ods are appropriate, errors introduced at the QM/MM 
boundary typically do not manifest themselves through- 
out the system. 

In metallic systems, however, the situation is quite dif- 
ferent. Bonds are not localized or associated with a dis- 
tinct pair of atoms. The embedded-atom picture^, IT(i| 
provides a more apt description of the situation. In 
the embedded atom picture of a simple metallic system, 
the density of the system is approximately the sum of 
charge densities of isolated atoms, and the energetic con- 
tribution of an individual atom to the system energy 
is approximately the embedding energy of the atom in 
a ho mog eneous electron gas. This picture, in various 
has been used to great effect to create 
classical pair functionals for metals. The success of these 
potentials in capturing the energetics of simple metallic 
systems, combined with their foundation on density func- 
tional theory arguments, make them ideal candidates for 
evaluating E 2 [II] in the present formalism. 

The notions of the embedded-atom method can be ex- 
tended to describe the energetics of a metallic region (re- 
gion /) within another metallic region (region II); region 
/ is embedded within region II. The exact nature of the 
embedding can be formally written in the manner of Eq. 
||TJ with density functional theory arguments. We first 
decide which ions will be associated with region I and 
which will be in region II, and we will denote those sets 
of nuclear coordinates by R 7 and R 11 . We denote the 
set of all nuclei with R tot = R 1 U R /7 . According to 
the Hohenberg-Kohn theorem, the total system energy, 
within the Born-Oppcnhcimer approximation, is given by 
minimizing a functional of the total charge density: 

E[R tot ] = mm E BFT [p tot ; R tot ] (3) 
p 

In order to be explicit, by Euft[p', {R-}] we mean: 
E DFT [p- R] = T s [p] + E u [p] + E xc [p] 

+ E/p(r)W-R ! )dr + E^ (4) 

i i<j 1 J 1 

where T s is the non-interacting kinetic energy functional, 
En is the Hartree energy, E KC is the exchange-correlation 
energy, and Vp Sp is the ionic pseudopotential. Thus 
Edft represents the combined electronic and ion-ion 
(Madelung) energy. 

If p tot is partitioned into two sub-densities, p 1 and p 11 , 
such that p tot = p 1 + p 11 , then the -Edft can be parti- 
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tioned: 

EuFT[p tot ', R tot ] = -EdftIp 7 ; R 7 ] 

+£dft [/' ; R /J ] + E int [p 1 , p 77 ; R 7 , R 77 ] (5) 

where i? mt is defined as in Eq. (J2J. By varying the to- 
tal energy with respect to p 7 , we find that the potential 
felt by p 1 is equal to the sum of the potential from re- 
gion I alone, plus an embedding potential V ela b(r) that 
completely represents the effect of region II upon region 
I: 

^ DFT [p tot ;R tot ] (5£ DFT [p 7 ;R 7 ] l , , 
J? = 6? + y ° mb(r) 

V emh (r) = ^V'Q*^ (6) 

By using different approximations for the terms in Eq. 
©, different coupled methods are obtained. Wesolowski 
and Warshel|l3j. building on the formalism of Cortona 
[l4| . used this partitioning to describe an efficient DFT 
method. In their scheme E[I] and E[II] are treated with 
Kohn-Sham DFT, but E int is evaluated with "orbital- 
free" density functional theory (OF-DFT), i.e. pure 
density functional theory in which the Kohn-Sham or- 
bitals are not used and the non-interacting kinetic en- 
ergy is app roximated with an explicit functional of the 

density|IalIlEllIllIIi- This allows E l T ] and E i n ] to 
be alternately minimized in the embedding potential of 

the other. Govind et al.0 utilized the partitioning of Eq. 

(J3J) to obtain a quantum chemistry(QC)/DFT coupled 

method. There Ei[I] was calculated with QC, E 2 [II] 

with DFT, and again E int was based on OF-DFT. They 

used this method to explore the electronic structure of 

molecules adsorbed on metal surfaces. Recently Kliincr 

et al. [2fJ have extended this formalism to treat adsorbed 

molecules in their excited state. 



III. FORMALISM 

The present method follows in the same vein as the 
last few examples, to achieve a DFT/classical potential 
coupling. The general idea of the present methods is 
as follows. Ei[I] is to be calculated with DFT. E 2 [II] 
is calculated via a classical potential. A choice can be 
made for the calculation of E mt , which results in distinct 
coupling methods, which we examine in detail below. 

A. Classical interaction energy 

E mt can be calculated using the classical potential: 

E int [I, II] = E cl [I + II] - E cl [I] - E cl [II] (7) 

Although this interaction energy is intended to represent 
the same DFT interaction energy that appears in Eq. 
0, it is not contradictory to use the classical potential 



to evaluate it, since the classical potential energy, evalu- 
ated for a given ionic configuration R, can be viewed as 
an approximation to the DFT functional that has been 
minimized with respect to the density; that is: 

E c \ [R] ~ min -Edft [p, R] ■ (8) 
p 

This choice of interaction energy results in a total energy 
of: 

E[R tot ] = £ cl [R tot ] - £ cl [R 7 ] + min£ DFT [p 7 , R 7 ] (9) 

p 1 

In this scheme, the forces on all atoms in region II are 
identical to forces on corresponding atoms if the classi- 
cal potential were used for the entire system; i.e. the 
DFT region has no effect on these atoms. If the cut- 
off length of classical potential is r c , then atoms that lie 
within region / and are farther than r c from the bound- 
ary will experience a force entirely from £<dft [I] j these 
atoms feel a force no different than corresponding atoms 
in a DFT calculation of region I. The force on atoms 
in region I that are within r c of the boundary do not 
come entirely from £"dft [I] , but also have contributions 
from E c \ [I + II] — E c \ [I] . These contributions should in 
principle be corrections to the surface forces experienced 
by these atoms from Edft[I]- Classical potentials have 
been developed to mimic the energetics, forces, and ge- 
ometries obtained from DFT calculations of various con- 
figurations, including surfaces 21]; such potentials should 
be particularly apt for the present coupling scheme. 

The implementation of this method demands nothing 
beyond what is required for a DFT calculation and a clas- 
sical potential calculation. It should be noted, however, 
that the DFT calculation, -Edft [I] , is a non-periodic cal- 
culation, and if OF-DFT is to be used, special consider- 
ations may need to be made for the calculation of non- 
periodic systems |23- 



B. Quantum interaction energy 

Alternatively E lnt can be calculated more accurately 
with a quantum mechanical method. Although we only 
represent region II by the coordinates of the ions of re- 
gion 77 atoms and calculate the energetics with a clas- 
sical potential, there is an implicit charge density p 77 
associated with £? c i[R 77 ] via Eq. JHJ. Because of this, 
we can consider a more sophisticated coupling scheme 
where the interaction energy is based on density func- 
tional theory. However, in order to compute the interac- 
tion energy via DFT when all we know about region 7T 
is an approximation of its charge density, the traditional 
Kohn-Sham scheme of DFT is not suitable. In the Kohn- 
Sham scheme, we start with a potential and obtain the 
density and energy of electrons in this potential. Instead, 
we need a means of calculating the energy of a system 
of electrons given their density. OF-DFT allows us to do 
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this. Thus we can write down the interaction energy in 
terms of OF-DFT energy functionals: 



E mt [I, II] = E OF [I + II] - E OF [I] - E OF [II] (10) 



At first glance this seems like a useless scheme, be- 
cause if DFT is used to calculate E mt [I, II], we may as 
well use DFT to calculate E[I + II], and thus no com- 
putational expense is saved with the coupled method. 
But because of the nature of many of the useful OF- 
DFT functionals, this turns out not to be the case. If 
E int [I, II] is calculated with OF-DFT, for typical approx- 
imate kinetic energy functionals the computation in Eq. 
(|10|l will require a computation time that is on the order 
of the computation time required to compute Eof [I] , the 
small subsystem, rather than the time required to com- 
pute Eof [I + II] . This is because significant cancellation 
is implicit in Eof [I + II] — Eof [II] ■ 

The existing approximate kinetic energy functionals 
differ in accuracy and computational efficiency. More- 
over, different choices of functional can be made for the 
evaluation of E[I] and E lnt , which further increases the 
number of possible coupling methods. This possibility is 
important because the degree to which the computation 
of E lnt can be made efficient depends on the choice of ki- 
netic energy functional and the functionals that will most 
efficiently treat E lnt are not necessarily accurate enough 
to treat the interactions within E[I]. 

Regardless of the choice of kinetic energy functional, 
the evaluation of E lnt [I, II] within this coupling scheme 
requires knowing the electronic density of region 77, 
p 77 (r). In the present method, /9 77 (r) is approximated 
as the sum of atomic charge densities p at (r) centered at 
the region II nuclei: 



nates within the method: 



P H (r) 



£p at (r-R 77 ) 



(11) 



£;[R tot ] = E cl [TL n ] +min 



£ F[p tot ;R tot ] 



Eqf[p ; R ] — Eqf[p ; R ] + £ DFT [p 1 ; R ] 



(12) 



The last term, -Edft [p 1 \ R 7 ] , is written as such (and not 
as Eof[I]) to emphasize that we could choose to com- 
pute it either with a Kohn-Sham-type calculation or with 
OF-DFT, but utilizing a more accurate kinetic energy 
functional than the other OF-DFT terms. This would 
allow for three distinct levels of accuracy in the calcu- 
lation: Kohn-Sham accuracy within region /, OF-DFT 
accuracy for the coupling between regions I and II, and 
the accuracy of the classical potential in region II. In 
this case, p 1 would consist of a set of Kohn-Sham or- 
bitals, p 1 (r) = |?/;i(r)| 2 , and we would minimize over 
the ipf. 



E[R tot ] = E cl [TL u ] 



mm 



£ F[p tot ;R tot: 



EofIp 11 -^ 11 ] - £ F[/;R I ]+£Ksh£> i ;R I ] 



P tot = £.N 

* — J i 

p 1 = E i^i 



(13) 

(14) 
(15) 



However, this interesting possibility is not explored 
presently; instead we use the same type of OF-DFT cal- 
culation for the last four terms of Eq. i|12|) . It should 
be noted that in this case the last two terms cancel, and 
then the total energy is given by: 



mm 

p 1 



E[H tot ] = E cl [TL H ] + 

£ OF [p tot ;R tot ]-£oF[p I7 ;R /z ] 



(16) 



This approximation is supported by the embedded-atom 
picture of simple metallic systems. In principle, p at (r) 
could be a non-spherically symmetric density. For ex- 
ample, if the arrangement of the region II atoms is al- 
ways close to the bulk lattice arrangement, then a non- 
spherically-symmetric charge density that reproduces the 
bulk charge density when periodically tiled may be more 
appropriate. However in this article p at (r) is always 
taken to be spherically symmetric. 

The density in region II is never explicitly represented 
in the calculation, but is given a precise form via Eq. 
(|ll|l . Thus region II is entirely described by the ionic 
coordinates R 77 , and p 77 , the form of which is needed to 
evaluate E lnt , is implicitly determined by R 77 . 

The second coupling method is summarized by the ex- 
pression for the energy as a function of nuclear coordi- 



C. Orbital-free DFT and approximate kinetic 
energy functionals 

Orbital-free DFT is a necessary part of the second cou- 
pling method, because the electronic structure of region 
II is represented only in terms of its density via Eq. 
(|ll(l : thus in order to utilize that information, E lnt must 
be based only on the charge density and the ionic coor- 
dinates. Here we describe some key ideas of OF-DFT. 

Hohenberg and KohnQ showed that the ground state 
energy of a system of electrons moving in an external po- 
tential is given by minimizing a density functional. Kohn 
and ShamQ wrote a useful partitioning of this energy 
functional: 

E[p] = T s [p] + E K [p] + E xc [p] + [ V e -i(r)p(r)dr (17) 
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where T s is the non-interacting kinetic energy functional, 
En is the Hartree energy, E xc is the exchange-correlation 
energy, and V e -i is the ionic potential. By introducing 
a set of fictitious non- interacting particles, we can ob- 
tain a set of single-particle equations, the Kohn-Sham 
equations, that allow for the evaluation of E[p] with an 
approximate E xc . The Kohn-Sham method results in an 
exact evaluation of T s [po] for the density po that mini- 
mizes E[p], but the method does not provide a means of 
evaluating T s [p] for an arbitrary density p. 

The Kohn-Sham partitioning of the energy, Eq. i|17[l. 
has turned out to be useful beyond the Kohn-Sham 
method. Because a number of limits of the exact non- 
interacting kinetic energy functional T s [p] are known |23j. 
T s [p] has been approximated by explicit density function- 
als constructed to satisfy one or more of these known 
limits. The orbital-free DFT methods are based on mini- 
mizing E[p] with T s replaced with an approximate kinetic 
energy functional. 

OF-DFT methods are typically more computationally 
efficient than the Kohn-Sham method. If the approxi- 
mate T s can be evaluated with an amount of computation 
that scales linearly with the system size, usually denoted 
by the total number of atoms N, then minimizing Eop[p] 
will require an amount of computation linear in the sys- 
tem size (O(iV) method). Since within OF-DFT all terms 
of the energy are explicit functionals of the density, there 
is no need for fictitious orbitals, and the density p(r) is 
the only represented variable. Thus, there is no need 
to solve the single-particle Scrhodinger equations for the 
fictitious particles while maintaining their orthogonality, 
operations which typically require most of the computa- 
tional effort in the Kohn-Sham approach and scale as a 
high power of the system size (0(iV 3 ) or higher). More- 
over, with the density p(r) as the only quantity of interest 
in the system, the OF-DFT methods use less memory 
than the Kohn-Sham method, since the latter requires 
the storage and update of a number of fictitious orbitals 
proportional to the system size, each of which consumes 
twice the storage (as complex quantities) needed for the 
density alone. 

In addition to computational advantages, unlike the 
Kohn-Sham method, the total energy functional Eof[p] 
can be evaluated for a given p(r). This property makes 
OF-DFT a suitable candidate for computing E^lp 1 , p H ] 
in the second coupling method discussed in the previous 
subsection. 

The number of available approximate kinetic energy 
functionals is sizeable, and the choice of functional is 
made based on considerations of efficiency and the types 
of systems to be treated. Because the systems to be con- 
sidered are simple metals with free-electron-like charge 
densities, an important property that should be included 
in the approximate kinetic energy functional is the cor- 
rect linear response around uniform densities: 



where f is the Fourier transform, and XLind(^) is the 
Lindhard response function: 



T 



<5 2 T S 



6p(r)5p(r>) 



XLindW — J 



1-X 2 

4x 



hi 



1 + x 



1-x 



(19) 



with k F — (Sir 2 po) 1 / 3 and x = k/2k F . 

A significant number of efficient functionals have been 
developed that satisfy the linear response limit for a par- 
ticular chosen average 

densitymdlHEEl- Such 
functionals often consist of several terms that are local or 
localized functionals, such as the Thomas-Fermi energy 
and the von Weizsacker functionals, plus one or more 
convolution terms: 



Tx\p] = jf(p(r))KQ 



XLind(fc) 



(18) 



r-r'|)ff(p(r'))drdr' (20) 



By choosing the kernel K(r) properly, the approximate 
functional can be made to satisfy the correct linear re- 
sponse, Eq. I|18|) . around some chosen uniform density po- 
Numerical tests indicate that among the current available 
efficient kinetic energy functionals, the ones of this form 
are most suitable for simple metallic systems. 

However, kinetic energy functionals that contain a con- 
volution part with a long-ranged kernel make the efficient 
evaluation of E mt [I, II] more difficult; the consequences 
of this will be discussed in the following section. 



IV. IMPLEMENTATION OF COUPLING 
A. Classical interaction energy 

The calculation of the energetics and ionic forces 
within the first coupling scheme described above involve 
only DFT calculations and classical potential calcula- 
tions. However, if an ionic relaxation is to be done on the 
whole system, there are several possible techniques, the 
optimal choice depending on the system being relaxed. 

If the partitioning of the system into regions / and 77 is 
such that the time required to calculate £<dft [I] is com- 
parable with the computation time of E c \ [I + II] , then 
ionic relaxation of the total system may be done by using 
a gradient-based minimizer such as conjugate gradients 
methods or quasi-Newton methods like BFGS 24] . If, on 
the other hand, the system partitioning is such that the 
time required to evaluate £7dft[/] is considerably more 
than that required for the computation of E C \[I + II], 
as is often the case, then an alternate relaxation scheme 
may be more efficient. The total system can be relaxed 
by using a gradient-based minimizer on the region I sys- 
tem alone, while fully relaxing the region // ions between 
each ionic update of region I. Gradient-based minimiz- 
ers like BFGS are only effective if the gradients involved 
are indeed gradients of an underlying object function. It 
is not immediately apparent that such is the case in this 
alternate-relaxation scheme, but we can demonstrate it 
as follows. 



The energy calculated with the first coupling scheme, 
as a function of all ionic coordinates, is given in Eq. JUJl- 
A secondary function that only depends on the region / 
ionic positions can be defined as: 



E'lH 1 



min MR*' 



(21) 



The useful aspect of E 1 is that its gradient with respect 
to Rf can be easily evaluated: 



dE' 



dE[K to 
<9Rf 

dE[R to 
9R| 



^ dE[R tot ] 9Rj; 



9Rf 



(22) 



(23) 



where the second term on the right of Eq. (1221) vanishes 
because all derivatives are evaluated at the minimum of 
_E[R tot ] with respect to R /7 . This result is analogous 
to the Hellmann-Feynman theorem[2!|. The introduc- 
tion of the E' function allows for the following relaxation 
algorithm: 

• minimize E[TL tot ] with respect to R H while holding 
R 7 fixed. This only involves the classical potential, 

#ci[R tot ]. 

• Calculate min p / E^ftIp 1 ', R 7 ] and £?ci[R J ], an d the 
forces on the region / ions. Using Eq. l|23[l the 
gradient of E 1 is obtained. 

• Perform a step of a gradient-based minimization of 
E'. 

• Repeat until the system is relaxed. 

In this manner, the number of DFT calculations per- 
formed is greatly reduced, albeit at the expense of more 
classical potential calculations. 



B. Quantum interaction energy 




FIG. 1: An illustration of the partitioning of the system ac- 
cording to the coupling method with quantum interaction en- 
ergy. 



The second coupling method maintains efficiency due 
to the cancellation that occurs when Eq F is computed. 
Consider the computation of a local term of E lnt ; such 
as the exchange-correlation functional within the local 
density approximation^ (LDA): 



= J /*c(p tot )dr-^/ xc (//)dr- J W)dr 
[/ xc (p tot ) - MP 1 ) - /xc(/ J )] dr (24) 



where / xc (p) = P e xc(p) and we have used the fact that 
p 7/ (r) = p tot (r) for r Q 1 . Thus calculation of E^ is an 
integral over f2 7 and not the entire system, which demon- 
strates our criterion for efficiency. Any local functional 
of p will obviously be calculated efficiently in the same 
manner. We note that when the same kinetic energy 
functional is used for the interaction energy and Eof[I] 
(which is the case for the tests performed in this paper), 
the cancellation exhibited in Eq. Ijltifl occurs. In this 
case, it is wasteful to compute the interaction energy as 
in Eq. 124|) and then compute and add on i? xc [p J ], as it 
exactly cancels the second term of Eq. <|24[) . Instead, we 
compute directly the following quantity: 



Implementation details of the second coupling method 
require more elaboration. One important point is that p 1 
must be confined to lie within a finite volume Q 1 . This re- 
gion should have significant overlap with the region where 
p 11 lies, in order to provide coupling of the two regions. 
But if p 1 were not confined to a finite volume fl 1 , it could 
in principle spread throughout the combined system, and 
during the course of minimizing with respect to p 1 (Eq. 
(|12|l h we would essentially be performing a DFT calcula- 
tion of the whole system. On the other hand, fl 1 should 
be chosen large enough so that p 1 is not artificially con- 
fined. In the test systems we examined, we found that 
when increasing the size of a point is reached where 
the results (e.g. the shape of p 1 and the forces on the 
ions) change little. The confinement of p 1 within ft 1 is 
illustrated in Fig. 2] 



E 



int+I 



[/xc(p tot ) - U C (P U )] dr 



(25) 



Similar considerations apply to the other parts of the en- 
ergy which are simple functionals of the density. The 
only term that does not fall in this category is the inter- 
action kinetic energy T s lnt , when it involves more sophis- 
ticated functionals with convolution terms such as Eq. 
(|2U|I . For this case, we have developed an appropriate ef- 
ficient methodology, the derivation of which is contained 
in appendix lAl 

Particular attention must be paid to the non-local 
terms of E lnt . As usual, cancellation occurs between 
electron-electron, electron-ion, and ion-ion terms that 
eliminates long-ranged interactions. The Hartree inter- 



7 



action energy is given by: 



pint 



p tot (r)p tot (r') - p J (r)/(r') - p H {v)p n {v') 



where f^ c = df xc /dp. For the long-ranged Coulombic 
functionals, the derivative is given by: 



drdr' 



drdr' 



/ /(r^V^r-R^dr 



where 



U H at (r) 



p at (r') 



dr' 



(27) 



Similarly the electron-ion interaction energy E^- reduces 
to: 

E^= f /(r)^F psp (r-Rf)dr 

Jn 1 



5> at (r-Rf)^V psp (r-R[)dr 



(28) 



where Vp Sp (r) is the pseudopotential representing the ion, 
and we have used Eq. to express p 11 (r) as a sum of 
p at . Finally the ion- ion interaction energy is given by: 



El 1 



E |R 



ZiZj 



Rfl 



(29) 



The combination of all three Coulomb terms can be ex- 
pressed as: 



E" 



El-\ = 



E< 



^ ^ ^clcc 

-R" 



) 



R 



dr 

(30) 



where we have defined: 



KM = V£\r) + V psp (r), 



Z. L Zj 



IR 



y psp (r-Rf)p at (r-Rf)di(31) 



Both V2g C (r) and 0y(R) are short-ranged functions in 
which the 1 / R dependence of the constituent terms can- 
cel. 

Within the second coupling method: (1) we minimize 
the energy with respect to p 1 , and (2) we calculate the 
forces on all of the ions and update their position. In or- 
der to minimize the energy with respect to p 1 , the deriva- 
tive SE int /Sp 1 '(r) needs to be calculated for r G fl 1 . This 
derivative can be evaluated efficiently for the local func- 
tionals like £v r : 



SE 



hit 

xc 



5/(r) 



/xc(p t0t )-/xc(/) 



(32) 



<V(r) 



E ™\ = F etec( r - Rf) (33) 



Evaluating this combined contribution to 5E mt /Sp 1 is a 
(2gjnrple matter of evaluating T/, a ' c for region 77 ions lo- 
cated near the boundary with region I. And so the gra- 
dient of the total energy with respect to p 1 is: 



SE 5J5 DFT [/;R J ] SE 



inl 



ST. 



hit 



Spi(r) 



<5/(r) 



<5/(r) Sp T {r) 



^^ C (r-Rf) (34) 



The calculation of the ionic forces proceeds differ- 
ently for region / and region 77 ions. Calculation of 
the region / ionic forces is facilitated by the Hellmann- 
Feynman theorem|2j|. If we denote the part of the en- 
ergy (Eq. H12|l 1 that is minimized with respect to p 1 by 
Gj/jB/.R 11 ]: 

G[/; H 1 , K H ] = E OF [I + II] - E OF [H] 

- E OF [I] + E DFT [I], (35) 

then we have, for the second coupling scheme: 

E[R tot ] = E cl [R n ] + minG[//;R / ,R // ] (36) 
p 1 

and when forces on region I ions are computed, the ex- 
pression simplifies: 



dE [TL b 
<9R{ 



dG 
dRf 
dG 
9RT 



6G d PLn( V ) 

Sp^r) dR( 



dr 



(37) 
(38) 



where the last term in Eq. (|37|l vanishes because we have 
minimized G with respect to p 1 , and so SG/Sp 1 ]^ = 0. 
So the forces on the region / ions are determined solely 
by the terms of G that explicitly depend on R 1 ; these 
terms are the electron-ion energy and the ion-ion energy. 
Using Eq. I|3(J|I . the force on the ith region / ion is given 
by: 

r _ dE\R™] 



d 



o 



dRf 

E C ^[I] + Ek-i[I\ + E^[I, II] + II] 



9R 



T [E e ^[I] + £U[J]] + £ V ^ ( R ' " R f ) 

3 



(39) 



Thus it can be seen from Eq. (|39|l that forces on region 
/ ions are given by the sum of the electron-ion and ion- 
ion forces present in subsystem / alone, and short-ranged 
interactions with region 77 ions that are nearby region 7 
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The forces on the region II ions come mostly from 
the classical potential, but they have contributions from 
E mt [I, II] because p 11 is a function of Rj 7 - Since we have 
not minimized with respect to p 11 , there is no Hcllmann- 
Feynman simplification when calculating the forces on 
region II ions, and all terms in the interaction energy 
contribute. The force on the jth region II ion is given 
by: 



dE[R tot ) 
dE cl [R n ] 



dE int [LII] 
d^ 1 



(40) 



Local functional parts of E mt such as the XC energy will 
have a contribution to the force given by: 



BE. 



int 



9Rf 



Vp at (r 



Rf)[/i c (p tot )-/i c (/ 7 )]dr (41) 



with analogous expressions for other local contributions 
that may exist in the kinetic energy functional such as 
the Thomas-Fermi energy. These local force contribu- 
tions are only non-zero for region II ions with an atomic 
density that extends into Q, 1 . It is also worth noting that 
the integral in Eq. I|41|) need not be carried out over all of 
f2 7 , but only over the intersection of f2 7 with p at (r— R 77 ). 

The Coulomb contributions to the region II ionic 
forces are given by: 



8 



int 



V0ij(Rf -R 



n 1 



/(r)VKt(r- 



R J )dr 



(42) 



This contribution also is non-zero only for region II ions 
near the boundary of il 7 . 

If a more sophisticated kinetic energy functional with 
a convolution term like Eq. (|20[1 is used in E lnt , then 
such a term adds considerable complication to the calcu- 
lation of the forces on region II ions, but these contribu- 
tions nonetheless die off as we move farther from region 
/. Thus within the framework of this coupling scheme, 
the forces on region II ions take the intuitively satisfy- 
ing form of being equal to the force that arises from the 
classical potential, plus a correction force for ions near 
the boundary of SI 7 . 

If the partitioning of the system between parts I and II 
is such that region I requires much longer computation 
than region 77, the second coupling method, like the first, 
allows for an efficient algorithm for ionic relaxation. We 
define a different partitioning of the ions as follows: we 
denote by R 7 the set of region / ions plus all region 
II ions R 77 for which the interaction force di? mt /<9R 77 

is not negligible, and we denote by R 77 the rest of the 




FIG. 2: A cut-away view of the test system. White atoms 
belong to region /, and dark atoms belong to region II. For 
the coupling method with quantum interaction energy, the 
region Q. 1 is shown by the white cube. 



R 77 ions. The point is that the forces on the R 77 ions 
only depend on the classical potential (as seen from Eq. 
(I40JI ) , and also that p 1 does not depend on the R 77 ions 
(as seen from Eq. l|H4|l). Thus the same algorithm for 
relaxing the system in the first coupling scheme can be 
used, but with R 7 replaced with R 7 , and R 77 replaced 
with R 77 . That is, before each relaxation step of the 
R 7 , the R 77 are to be fully relaxed. 



V. TESTS 

In order to test the present coupling methods, we have 
focused on a simple coupled system that is readily ana- 
lyzed. The system consists of 10 x 10 x 10 cubic unit cells 
(4 atoms each) of crystalline fee aluminum. The inner- 
most 2x2x2 unit cells (32 atoms total) are considered 
to be region /, and all atoms outside are considered to be 
in region II. Region II, which is treated with the classi- 
cal potential, is treated as a periodic system in order to 
remove effects of surfaces from the simulation. So in fact 
the test system consists of an infinite array of 32-atom 
Al clusters treated quantum mechanically, embedded in 
an Al bulk treated by classical potentials. Obviously, if 
there is good coupling between region I and region II, 
the entire system should simply behave like pure bulk fee 
Al, making it easy to evaluate the quality of the coupling. 
This test system is illustrated in Fig. [21 

In all of the present tests, region I is treated with 
OF-DFT. However, the particular kinetic energy func- 
tional used differs among the tests. The Al ions 
are represented with the Goodwin-Needs-Heine local 



pseudopotential 26] . For all tests, the classical potential 
used is the "glue" potential of Ercolessi and Adams [Tlj. 
which has the embedded-atom method (EAM) form: 



E[R t ]=J2F 



P 



EAM/ 



|R«I) 



5$>(|R«|)(43) 



The EAM potential has been scaled both in r and in 
energy: 



F[p] -> aF[p], 
3 eam (r) ^ p EAM (/3R), 

0(R) -> a0(/3R) 



(44) 



with a and /3 chosen so that the potential yields pre- 
cisely the same lattice constant and bulk modulus of fee 
Al simulated with OF-DFT employing the particular ki- 
netic energy functional used in that test. This is in accord 
with the philosophy that the coupled simulation should 
behave as if the accurate method were used for the entire 
system. But this procedure is also done so that a "fair" 
assessment of the coupling itself can be made; we wish 
to examine errors in the present coupling methods them- 
selves and the approximations involved in them, but not 
the errors coming from a trivial incompatibility between 
energy methods. To make the classical potential even 
more compatible with the OF-DFT method, we could 
re-determine the form of the classical potential using the 
method that Ercolessi and Adams originally used to de- 
velop their potential|2l||: perform a large number of ref- 
erence energetic calculations of Al using OF-DFT, and 
find the EAM-type potential that best reproduces these 
results. This would be a rather involved procedure, so 
we have chosen to simply scale the potential. 



A. Test of classical interaction energy method 




FIG. 3: Test of coupling method with classical interaction en- 
ergy, (a) The forces on the region / atoms when the atomic 
positions are at the perfect lattice positions. The force fac- 
tors are scaled so that a vector of length one lattice constant 
corresponds to 1 eV/A. (b) The relaxed region I and 2 atomic 
positions shown in white and black, respectively. The perfect 
lattice sites are drawn as gray spheres of a slightly smaller 
radius. 



In the first coupling method (Eq. @), the energetics 
of region I was treated with OF-DFT, and the kinetic 
energy functional used was one developed by Wang et 
al.yjj, with a density-dependent kernel, and parameters 

{a, (3, 7, p*} = {5/6 + V5/6JS/6 - Vb/6, 2.7, 0.183 A" 3 } 
(in the notation of Ref. [ljj ). This functional has six 
convolution terms of the form of Eq. (|20|l . The clas- 
sical potential was scaled to match the lattice constant 
(4.033 A) and bulk modulus (55.7 GPa) of fee Al ob- 
tained by this kinetic energy functional. 

The system was initially arranged in the perfect fee lat- 
tice configuration The forces on the region atoms were 
identically zero, since they come entirely from E cl [I +11], 
which is at a minimum in the initial configuration. How- 
ever, the forces on the region I atoms are not zero, as the 
OF-DFT and EAM forces do not perfectly cancel each 
other. The average magnitude of the force difference per 
atom between the OF-DFT and EAM calculations of re- 
gion / was 0.34 eV/A. These initial forces on the region 



/ atoms are shown in Fig. OJa), with the drawn force 
vectors scaled so that a force of 1 eV/Awould extend one 
lattice constant. Then the coupled system was relaxed. 
After relaxation, it was found that the atomic positions 
deviated from the correct fee lattice positions by an av- 
erage of 0.004 Aper atom. The average deviation of just 
the region I atoms was 0.076 Aper atom. The atomic de- 
viation is shown in Fig. E^b) , in which the relaxed atomic 
positions for the region / and region II atoms are drawn 
as white and black balls, respectively, and the correct lat- 
tice positions are drawn as gray balls of a slightly smaller 
radius. Note that only the four (100) layers that include 
region / are shown. From this diagram it can be seen that 
in general the relaxed atomic positions deviate from the 
the perfect lattice positions by bulging out from region 
/ slightly, with the deviation decreasing with increasing 
distance from region /. 
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B. Test of quantum interaction energy method 

The second coupling method was applied to the same 
simple test system. The kinetic energy functional em- 
ployed was again from Ref. (19]), but in this case it was 
a functional with a density-independent kernel, with pa- 
rameters {a, (3} — {5/6 ± v5/6}. A different functional 
was chosen for this test because of its simpler form: it 
contains only one convolution term of the form of Eq. 
(|20|l . while the functional used in the test of the first 
coupling method had six. This makes the evaluation of 
the kinetic interaction energy, T s mt , simpler. Further- 
more, this functional performs well for structures that do 
not deviate much from the bulk system. This functional 
was found to be inapplicable to the test of the first cou- 
pling method, because in that approach the calculation 
of Eof[I] amounts to an isolated cluster; in that case, 
there is no embedding potential from region // and the 
minimization with respect to p 1 does not converge. For 
bulk fee Al, however, this simpler functional employed to 
test the quantum interaction energy method, performs 
quite well producing an equilibrium lattice constant of 
4.035 Aand a bulk modulus of 71.9 GPa. 

Another aspect of the second coupling method is the 
choice of atomic density functions p at representing p 11 
through Eq. (|11|) . Two different choices of p at were 
tried. One choice, p at >s as j W as the valence density from 
a Kohn-Sham calculation of an isolated Al atom, rep- 
resented with the same Pseudopotential|26| . The other 
choice, p at < cr y st was again a spherically symmetric charge 
density chosen such that the charge density that results 
from periodically superposing it on an fee lattice most 
closely matches the charge density coming from an OF- 
DFT calculation of bulk fee Al. The desired spherically 



symmetric charge density p 



,at,cryst 



)(r) minimizes: 



f ^fec 
JO. 



p at < cryst )(r) -p fcc (r)] dr 



(45) 



where S (r) is an infinite fee lattice of delta functions, * 
denotes convolution, and is one unit cell; p (r) is the 
valence charge density of the fee crystal. In reciprocal 
space, this becomes: 



Q 



S Q P- 



~at,cryst 



(Q) 



Pq 



(46) 



where Sq is the structure factor. This is minimized by 
requiring: 



~at,cryst(Q) 



</5f cc 



Q IQ 



(Sq)c 



(47) 



where (• • • )q denotes an averaging over reciprocal lattice 
vectors Q of length Q. This p at ' cryst (Q) was then used 
to construct a radial charge density p at > cr y st (r) that was 
commensurate with p at,cryst (Q) at the values of Q where 
the latter was defined. The two charge density choices 
are shown in Fig. 



•< 




FIG. 4: The two atomic densities used to represent p 1 



The second coupling method was tested on the same 
system used to test the first coupling method. With the 
second method, however, we must choose the form of the 
atomic density representing p 11 , and the extent of the 
region fl 1 . With respect to the choice of f2 7 , we have 
found the following general behavior: if fl 1 is chosen to 
be too small, then after minimization with respect to p 1 , 
there is an excess buildup of charge near the boundary of 
ft 1 . This in turn results in a net attraction of the region 
I ions toward the boundary of £1 7 . This is remedied by 
an increase in the size of il 1 . When is increased fur- 
ther still, the results (the ionic forces, and p 1 ) are found 
to change only very slightly. We note that regardless of 
the size of Q, 1 , the total density is always found to be 
continuous at the boundary due to the high energy that 
the kinetic energy functional assigns to a discontinuity 
in the density. In Fig. (JSJ we have plotted the total 
charge density after minimizing with respect to p 1 , with 
p 11 given by (a) a superposition of p 3 *'®*, and (b) a su- 
perposition of p at < cr y st . The particular slice of the charge 
density is a (100) plane that passes through one of the 
central atomic planes of the region / cluster. In (c) and 
(d) we have plotted the difference between these coupled 
charge densities and the density of this system when com- 
puted entirely with OF-DFT. In general, using p at . cr y st 
results in a more accurate total charge density. From (c) 
and (d), it is clear that the superposition of p at > cr y st re- 
produces the pure OF-DFT crystal charge density better 
than p at >s as both for points r well within Q, 1 , as well as 
at the boundary of Vl 1 . 

It turns out, however, that the forces on the ions, for 
both choices of p at , are comparable. Also comparable is 
the amount of deviation from the perfect lattice positions 
upon atomic relaxation for both choices. The exact num- 
bers for these quantities for both coupling methods and 
the two choices of p at are summarized in Table [I] 
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FIG. 5: Results for the quantum interaction energy method, (a) and (b) are contour plots of p tot with p 11 given by superpositions 
of p at ' gas and p at - cr y st ) respectively. The boundary of Q 1 is shown with a dashed line, and the positions of the region 7 atoms 
lying in this plane are indicated by (+), and region 77 atoms by (x). (c) and (d) show the difference between these two densities 
and the "correct" density coming from a purely OF-DFT calculation of the whole system. 



TABLE I: Summary of the performance of the two coupling 
methods, and the two choices for p at in the quantum interac- 
tion energy method. F„ mx , F^f ax are the maximum forces on 
region I atoms and region 77 atoms, and dmaxi t the max- 
imum displacements from the perfect lattice positions upon 
relaxation for region 7 and region 77 atoms. 



interaction 


F 1 

■*■ max 


F 11 

± max 


^max 


d n 

^max 


energy 


(eV/A) 


(eV/A) 


(A) 


(A) 


classical 


0.45 





0.12 


0.048 


quantum, p at > gas 


0.12 


0.27 


0.12 


0.15 


quantum, p**^* 


0.094 


0.37 


0.217 


0.265 



VI. CONCLUSIONS 

The coupling of classical and quantum simulation in 
simple metals involves a set of challenges quite different 
than those for the coupling of covalently-bonded materi- 
als and molecules, and hence requires a different set of 
approaches. We have presented here two methods for 
combining classical and quantum mechanical simulation 
of simple metals. Both are based on a similar partition- 
ing of the energy of the system, but they differ in how the 
energy of interaction between the classical and quantum 
mechanical parts of the system are treated. We have pre- 
sented numerical implementations that allow both cou- 
pling methods to be efficient. 

Within the first coupling method, in which the inter- 
action energy is determined from the classical potential, 
forces in the classical region are fully determined by the 
classical potential. Forces in the quantum region are de- 
termined by both classical and quantum energetics, the 
quantum energetics dominating well within the quantum 
region. A major practical advantage of this approach is 
that, if region / contains many different atomic species 
while region 77 contains only one atom type, there is no 
need for a classical potential for each species and their in- 
teractions; if the various species of atoms are well within 
region /, then the potential representing them does not 
matter at all as interactions in this region are treated 



purely with quantum mechanics. 

Within the second coupling method, in which the in- 
teraction energy is determined via OF-DFT, forces in the 
classical region are determined mostly by the classical po- 
tential, with quantum contributions to atomic forces near 
the boundary of the regions. Forces in the quantum re- 
gion are determined fully by quantum energetics. Within 
the quantum region, the charge density accurately repro- 
duces the correct charge density, and smoothly joins with 
the implicit density (given by a sum of atomic densities) 
of the classical region. 

Test results indicate that the second coupling method 
yields more accurate forces on the atoms in the quantum 
region than the first method, but that the first method 
yields more accurate forces for the atoms in the classical 
region. This may be due, to some extent, to the less- 
accurate OF-DFT method used in the test of the second 
coupling method. The first coupling method also yielded 
a better relaxed structure, probably due to its better 
treatment of forces on the classical atoms. However, un- 
like the first method, the second coupling method results 
in a more accurate charge density within the quantum 
mechanical region, allowing for an accurate treatment of 
physical problems such as the introduction of impurities, 
where the background density is important. We also find 
that a superposition of atomic charge densities can repro- 
duce the actual charge density well for a simple metallic 
system, given an appropriate choice for the atomic charge 
density; this allows for a smooth density transition at the 
boundary between regions. 

Clearly, an important issue affecting the coupling qual- 
ity for both methods is the agreement between forces 
from the DFT methods; within both methods there are 
atoms whose forces are determined by a combination of 
quantum and classical energetics, and the more closely 
the two energetics agree, the better the coupling will be. 
An improvement in the quality of the coupling might 
be obtained if the classical potential employed in region 
II is optimized to closely reproduce the DFT energetics; 
this is also in accord with the multiscale philosophy that 
a coupled simulation should act as if the most accurate 
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method were used to simulate the entire system. 



APPENDIX A: EVALUATING THE 
INTERACTION ENERGY FOR COMPLEX 
KINETIC ENERGY FUNCTIONALS 

We describe here how the interaction energy can be ef- 
ficiently calculated when the approximate kinetic energy 
functional used is of a more complicated form, containing 
a convolution term of the form of Eq. (|20|l . That is, we 
will describe a method for evaluating: 

2kV,/>"] = / /i 2 (r)i^r-r')5i2(r')drdr' 
-J h{v)K{v~v') 9l {T')Avdv' 
- J / 2 (r)X(r - r'). 92 (r')drdr', (Al) 

where we have defined /i(r) = /( j o / (r)), /i2(r) = 
f(p I (r) + p 11 ^)), and so on. Then we define two new 
functions, 



F(r) ee / 12 (r) - / a (r), 
G(r) ee g 12 (r) - g 2 (r) 



(A2) 



Note that F(r) and G(r) are zero for points r £ Vl 1 . 
Using F and G we can re-express Eq. (JAlf) as: 

T^= [ F(r)(X*G)(r)dr- f h(v)(K * ffl )(r)dr 

+ f F(r)(K * 52 )(r)dr + / G(r)(K * / 2 )(r)dr 
Jn 1 Jo. 1 



(A3) 



where we have now defined: 



(K * G) (r) ee /" K(r - r')G(r')dr', etc. (A4) 

We point out that if this interaction energy is being cal- 
culated in a coupled simulation in which the energetics 
of region / are calculated using the same kinetic energy 
(i.e. E[I] and E lnt [I,II] being treated at the same level 
of theory), then the final term of Eq. (|A3|) is equal to 
and will cancel with the corresponding term in E[I]. 

So with Eq. (|A3(I we have expressed T^ 1 ' purely in 
terms of intergrals over fl 1 ; the problem is now reduced 
to efficiently calculating the functions (K * j^Xr) and 
(K * g 2 ) (r) for points r within ft 1 . A straightforward 
integration for each point r G is not an option, be- 
cause K(r) is typically long-ranged, and thus determin- 
ing (K * / 2 )(r) at one single point r would require an 
integration over the volume of the whole coupled system, 
which would be highly inefficient. We now describe a 
method for determining (K * f 2 )(r), and (K * g 2 )(r) can 
be determined with precisely the same method. 



In earlier work_22] we have developed a method for ef- 
ficiently evaluating convolutions like Eq. I|A4(I when the 
convolution kernel K(r) is of the particular form typi- 
cally encountered in kinetic energy functionals involving 
convolution terms |15l Il6t Il7l ll8L Il9|| . We will invoke this 
method to determine (K*f2). In this method, the kernel 
is fit in reciprocal space with the following form: 



K(k) ~ Y^kiik), 

i 

P k 2 



k 2 + Qi 



(A5) 



where Pi, Qi are complex fitting parameters. The ker- 
nels encountered in many kinetic energy functionals are 
well-fit with this form, with 4 terms. The kernel in real 
space can be expressed as the sum of the inverse Fourier 
transform of each term of Eq. (|A5|) : 



K(r) ~ 



Ki(r) ee Pi5(T)-PiQi 



(A6) 



Thus {K * f 2 ) can be written as the sum of separate 
convolutions: 

(K*f 2 )(r) = ]T(/^*/ 2 )(r), (A7) 

i 

(Ki*h)(T) ee y^(r-r')/2(r')dr' (A8) 

Because in reciprocal space the [Ki * / 2 ) satisfy: 

[k 2 + Qi] (Jf7T/ 2 )(k) = i^ 2 / 2 (k), 
in real space they satisfy: 

[V 2 -Q a ] (JWa)(r) =^V 2 / 2 (r) 



(A9) 



(A10) 

i.e. they are solutions to (complex) Helmholtz equa- 
tions which can be solved with conjugate-gradient-based 
methods |27j: such methods are efficient and only involve 
operations within fl 1 . The solutions to Eqs. I|A10|) are 
only well-defined when boundary conditions for (Ki * 
/2)(r) are supplied. 

We propose the use of Dirichlet boundary conditions. 
The value of (Ki * /2)(r) for points r on the boundary 
of Q 1 can be found by evaluating the convolutions, Eqs. 
(|A8|) . Because of the regular nature of p 11 (r), being the 
sum of atomic densities, an efficient method for evaluat- 
ing these convolutions exists. The form of the convolu- 
tion that needs to be evaluated is: 



(#i*/2)(r) 



J Ki(r — v')f ^> at (r'-Rf)j dr' 



( AU ) 

If f(p) were a linear function, then this would reduce to 
a sum of pair functions. For many kinetic energy func- 
tionals, f(p) is not linear, but equal to some power of 
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1.5 



0.5 - 



f(p)=P" 

2nd-order f-expansion 

1 st-order h-expansion 





FIG. 7: New coordinates R c and R rc i for evaluating the three- 

(2) 

center integrals of L\ . 



FIG. 6: A demonstration of the Taylor expansion of Eq. 
1A12II compared to a direct Taylor expansion of f(p) = p 1 ' 5 
about po = 1. 



P : f(p) — P a ■ This leads us to consider a Taylor expan- 
sion of f(p) about some average density p Q . This Taylor 
expansion suffers in places where p 11 (r) is near 0, which 
occurs, for instance, in the center of fl 1 . An expansion 
that is much more accurate down to small values of p is 
obtained if we Taylor expand the function h(p) = f{p)/p 
and express f(p n ) in terms of this expansion: 



h(p ) + h'(po) 



5> at (r'-Rf)- 



Po 



A12) 



In Fig. El we illustrate the effectiveness of l|A12j) com- 
pared to expanding f(p) directly when f(p) — p and 
Po = 1. Although the h-expansion is taken only to first 
order, while the f-expansion is taken to second order, the 
h-expansion is seen to be more accurate at small p. 

Upon substitution of the expansion of Eq. I|A12(I in 
the convolution, Eq. i|All() . we find: 



if } (R,R') ee 



[h(p )-p h'(p )]^2L ( i 1 \r-Jiy) 

3 

+ ^o)$>f ) (r-Rf,r-Rf), 

J X 4 (r> at (R-r')dr', 

J Ki(r')p at (R - rV*(R' - r')dr' 

(A13) 



is the convolution of an atomic density with Ki(r). 

( 2) 

L\ (R, R') is the convolution of the product of two 
atomic densities with K(r), and consequently vanishes 



when the two atomic densities do not overlap. The inte- 
grand is non-zero only where the overlap occurs. It thus 

(2) 

makes sense to express L\ in terms of new coordinates, 
illustrated in Fig. [7| 

Mi(R c ,R re i) = Lf >(Re - ±R rcl , R c + |R rel ) (A14) 

In the case of a spherically symmetric p at (r), L^^R) 
is a function only of |R|, and Mj(R c ,R rc i) will depend 
only on |R C |, |R re i|, and R c • R rc i. We now argue that 
the dependence on R c • R rc i is weak. We can write the 
expression for Mi as: 

Mj (Rc, R re i) = J KiiRc - r)p at (r - iR rcl ) 

xp at (r+iR rcl )dr (A15) 
If we expand K about R c , we find: 

M 4 (R c ,R rcl ) = J [K i (\K c \)-r-VK i (R c ) + ---} 

x p at (r - |R re i)p at (r + |R re i)dr (A16) 

The integrals over the odd powers of r in the expansion 
of Ki vanish by symmetry. Thus if we truncate the series 
at first order, the first order term vanishes, leaving only 
the 0th order term: 

MiCRcRrei) ~ A"j(|R c |)P(|R rcl |), (A17) 
P(|R rcl |) = y"p at (r)p at (r-R rol )dr (A18) 

Truncating the expansion of Ki at the first order is rea- 
sonable, because Ki(r) oscillates around the Fermi wave- 
length of the system, a length scale close to that of p at . 
This approximate form, Eq. 1A17|I will behave quite 
badly at small i? c , because Ki diverges at the origin, 
and the radial averaging of this divergence that occurs 
in Eq. HA15fl is not reflected in Eq. IJA17II . Thus we 
replace Xi(|R c |) there with the convolution of Ki with a 
Gaussian of unit weight and a variance ro given roughly 
by the length-scale of the overlap regions of the atomic 
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densities (i.e. some fraction of the range of /? at (r)): 

M,(R c ,R rcl ) ~ iv'(|R c |)P(|R rcl |), 
K[(v) = (K^w)(r), 

w(t) ee Tr- 3 / 2 r- 3 e-( r/ro)2 (A19) 

Summarizing these results, we find that we can ap- 
proximately evaluate (Ki * f 2 ) as: 

(Ki * / 2 )(r) ~ [h(p ) - poti(p )] J2 - Rf ) 

j 

+ h'(po) £ ^'(|r-I(Rf +R»)\)P(\B? -Rf|) 

<j,k> 

(A20) 

where the summation over < j,k > indicates that we 
need only sum over pairs of region II atoms with over- 
lapping densities. The derivation of Eq. I|A20|) involved 
several approximations, and thus is not expected to be 
precise. We only propose that Eq. (|A20|I be used to 
generate the boundary conditions for Eqs. IjAlOjl that 
determine the (Ki * /2)(r) within region f2 7 , and we have 
found that the resulting (Ki * /2)(r) is more dependent 
on the source term than the boundary conditions. Nev- 
ertheless, because of the inaccuracies of Eq. (| A20I) . we 
define a new region, ft 1 , that contains and extends a bit 
beyond ft 1 , and we use Eq. I|A20|I to obtain the bound- 
ary conditions for points r that lie on the boundary of 
f2 7 , and we solve Eqs. (|A10|) for all r S , so that the 
resulting (Ki * /2)(r) are accurate for all red 1 . 

Thus we have all the pieces necessary to compute Tj?* 
In summary, we do this as follows: 



• Using Eq. I|A20|) . we can evaluate (Ki * / 2 )(r) and 
(Ki*g 2 )(r) for points r on the boundary of a region 
SV that is slightly larger than tt 1 . 

• Using those boundary conditions, the Helmholtz 
equations 1A10|I are solved, yielding (Ki * / 2 )(r) 
and (Ki * 32) (r) for all points rgtf . 

• Then (K * /2X1") and (K * g 2 ) (r) are constructed 
with Eq. 1A7(I . and we can evaluate T^ 1 * via Eq. 

Ther kernel interaction energy T^ 1 * also gives a small 
contribution to the forces on region II atoms near the 
1-2 boundary. By differentiating Eq. l|Aljl . we find: 

^77 = - J Vp(r-Rf)[# 3 (r)(ir*G)(r) 

+g' 12 (r)(K * F)(r) + F'(v)(K * g 2 )(v) 
+G'(r)(X*/ 2 )(r)]dr, (A21) 



(K*F)(r) = J K(r-r')F(r')dr', 

f[ 2 (r) ee r(p I (r) + p H (r)), 
F'(v) ee f 12 (v)-f(p u {v)), etc. 
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